#ifndef DYN_ODE_SIFS8_H
#define DYN_ODE_SIFS8_H

#include "dyn_ode_si.h"

/*! \brief Symplectic integrator with fixed stepsize of order 8.
 *
 * We implement the eighth-order
 * symplectic integrator that Yoshida found.
 * with his elegant method.
 *
 * This integrator only works with Hamiltonians of the
 * form
 *  \f[ H(\mathbf{q},\mathbf{p}) = T(\mathbf{p})+V(\mathbf{q}), \f]
 * i.e., separable Hamiltonians.
 */

class SIFS8 : public SI
{
public:
    /*! Same as base constructor. */
    SIFS8(HamiltonianSystem& func,
          colvec _initCond,
          double _end,
          double _initStepsize)
        : SI(func,_initCond,_end,_initStepsize), M(16), order(8){}

protected:
    /*! @name Coefficients of the expansion. Solution D of Yoshida (1990). */
    //@{
    // Auxiliary variables
    const double w1 = 0.102799849391985e0;
    const double w2 = -0.196061023297549e1;
    const double w3 = 0.193813913762276e1;
    const double w4 = -0.158240635368243e0;
    const double w5 = -0.144485223686048e1;
    const double w6 = 0.253693336566229e0;
    const double w7 = 0.914844246229740e0;
    const double w0 = 1.0-2.0*(w1+w2+w3+w4+w5+w6+w7);

    // Actual variables
    const double c1 = 0.5*w7;
    const double c2 = 0.5*(w7+w6);
    const double c3 = 0.5*(w6+w5);
    const double c4 = 0.5*(w5+w4);
    const double c5 = 0.5*(w4+w3);
    const double c6 = 0.5*(w3+w2);
    const double c7 = 0.5*(w2+w1);
    const double c8 = 0.5*(w1+w0);
    const double c9 = c8;
    const double c10 = c7;
    const double c11 = c6;
    const double c12 = c5;
    const double c13 = c4;
    const double c14 = c3;
    const double c15 = c2;
    const double c16 = c1;

    const double d1 = w7;
    const double d2 = w6;
    const double d3 = w5;
    const double d4 = w4;
    const double d5 = w3;
    const double d6 = w2;
    const double d7 = w1;
    const double d8 = w0;
    const double d9 = d7;
    const double d10 = d6;
    const double d11 = d5;
    const double d12 = d4;
    const double d13 = d3;
    const double d14 = d2;
    const double d15 = d1;
    const double d16 = 0.0;

    const colvec c = {c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16};
    const colvec d = {d1,d2,d3,d4,d5,d6,d7,d8,d9,d10,d11,d12,d13,d14,d15,d16};
    //@}
};

#endif // DYN_ODE_SIFS8_H
